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1  Introduction 


Improvised  explosive  devices  (IEDs)  oftentimes  hidden  along  roads  and  detonated  underneath 
trucks  or  tanks  pose  a  major  threat  to  soldiers  and  can  cause  significant  accelerative  injuries.  The 
explosion  beneath  the  vehicle  is  commonly  referred  to  as  an  underbody  blast  (UBB)  event.  It  has 
been  reported  that  8,159  IED  incidents  occurred  in  Afghanistan  in  2009,  which  was  an  increase 
from  3,867  and  2,677  which  were  reported  in  2008  and  2007,  respectively  [Whitlock,  2010].  While 
civilian  car  traffic  accidents  are  in  the  range  of  5  to  35  G  [Gabauer  and  Gabler,  2008],  an  underbody 
blast  may  cause  accelerations  in  the  lower  extremities  in  the  range  of  155  to  217  G  [Nilakantan  and 
Tabiei,  2009].  During  an  underbody  blast  event,  within  0.5  ms  of  the  initiation  of  the  explosion  a 
shock  wave  hits  the  bottom  plate  of  the  vehicle  and  causes  a  rapid  rise  in  pressure.  This  pressure 
results  in  local  acceleration  and  subsequent  deformations  of  the  plate  which  apply  significant  loads 
to  the  soldier’s  lower  extremities  [NATO  HFM-090,  Task  Group  25,  2007].  As  the  blast  wave  re¬ 
flects  under  the  vehicle,  a  pressure  force  acts  on  the  bottom  of  the  vehicle  and  causes  it  to  lift  off  the 
ground.  The  acceleration  of  the  vehicle  and  the  collisions  that  follow  are  also  a  significant  source 
of  injury,  especially  if  the  soldier  is  not  properly  restrained.  In  summary,  the  body  is  subjected 
to  both  the  local  effects  (caused  by  shock  and  high  rate  deformation  of  the  vehicle)  and  global 
effects  (vehicle  motion  over  time)  as  a  result  of  the  mine  detonation  process.  The  local  effects  are 
less  understood  and  lead  to  a  range  of  injuries  including  fractures  of  the  knee-thigh-hip  complex, 
calcaneus,  ankle,  midfoot,  tibia  and  fibula  shaft  and  malleoli,  as  well  as  knee  and  other  ligament 
tears  and  ruptures. 

Some  of  the  understanding  for  injuries  experienced  within  the  underbody  blast  environment 
can  be  obtained  from  the  civilian  vehicle  crash  environment.  For  example,  fractures  of  the  ankle 
and  lower  foot  which  are  predominant  within  the  military  environment,  have  already  been  studied 
to  some  degree  in  the  past.  Manning  et  al.  [1998]  report  that  calcaneal,  talar,  midfoot  and  various 
other  ankle  fractures  have  been  attributed  to  the  mechanism  of  axial  loading  through  the  plantar 
surface  of  the  foot.  While  these  injuries  might  be  expected  from  axial  loading,  Crandall  et  al. 
[2003]  demonstrated  that  even  malleolar  fractures  can  occur  from  pure  dynamic  axial  loading  of  the 
leg.  Begeman  and  Paravasthu  [1997]  also  conducted  experimental  impact  tests  where  unembalmed 
cadaveric  legs  were  subjected  to  a  uniaxial  plantar  surface  load  along  the  axis  of  the  tibia.  Pilon 
and  calcaneal  fractures  were  observed  with  an  average  tibial  axial  force  of  7590  N  at  failure.  Kuppa 
et  al.  [2001]  offers  a  useful  review  of  calcaneus,  talus,  ankle  and  midfoot  fracture  research,  most 
of  which  was  conducted  on  postmortem  human  subjects.  While  this  past  research  highlights  some 
common  features  with  military  blast  injuries,  it  is  important  to  note  that  the  rate  of  loading  to  the 
lower  extremities  may  be  10  -  100  times  greater  than  civilian  vehicular  crashes.  Thus,  there  is  a 
strong  possibility  that  different  injury  mechanisms  may  be  activated  during  these  events  and  new 
or  enhanced  injury  criteria  may  be  required. 

There  have  been  attempts  to  computationally  model  underbody  blast  effects  in  the  lower  ex¬ 
tremities.  Nilakantan  and  Tabiei  [2009]  attempted  to  assess  injury  by  measuring  tibia  compressive 
loads  following  an  underbody  blast  event  using  finite  element  analysis  of  a  HYBRID  III  human  sur¬ 
rogate.  They  make  clear  that  more  advanced  and  biofidelic  lower  extremity  finite  element  models 
are  needed.  Many  lower  extremity  studies  employ  a  combined  experimental  and  simulation-based 
approach.  For  example,  Manseau  and  Keown  [Manseau  and  Keown,  2005]  compared  a  physical 
lower  leg  model  and  a  computer  model  in  the  development  of  enhanced  injury  assessment  methods, 
with  a  focus  on  the  ankle  complex  and  its  connections  to  the  tibia.  Bandak  et  al.  [2001]  developed 
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an  anatomically  realistic  finite  element  model  of  the  lower  leg  to  study  the  effects  of  axial  impact 
loading  to  the  plantar  surface  of  the  foot.  The  model  included  explicitly  represented  bones,  cortical 
and  trabecular  bone  types,  and  various  soft  tissues.  The  plantar  surface  of  the  foot  was  impacted 
by  a  rubber-layered  plate  within  a  range  of  velocities  up  to  6.7  m/s.  The  model  was  compared  to  a 
set  of  experiments  in  which  human  cadaveric  lower  legs  were  impacted  by  a  swinging  pendulum 
at  various  initial  velocities  ranging  from  2.2  m/s  to  7.6  m/s  [Yoganandan  et  al.,  1995,  1996].  In  the 
experiments,  bone  fracture  was  observed  in  the  distal  tibia  and  calcaneus  at  velocities  at  and  above 
6.7  m/s.  The  simulations  by  Bandak  et  al.  [2001]  compared  reasonably  well  with  the  experimental 
data  of  Yoganandan  et  al.  [1995,  1996],  though  an  over-prediction  in  peak  load  values  at  the  proxi¬ 
mal  end  (i.e.  knee)  and  plantar  surface  of  the  foot  became  increasingly  pronounced  with  increasing 
velocity  (4.47m/s  and  above),  perhaps  because  failure  was  not  incorporated  into  the  finite  element 
model. 

In  this  paper  we  develop  a  numerical  model  of  the  lower  leg,  constructed  in  an  effort  to  un¬ 
derstand  the  dynamic  failure  of  the  lower  extremities  when  subjected  to  underbody  blast  loading. 
This  body  region  (from  the  knee  down)  is  the  area  most  often  injured  in  an  underbelly  blast  event 
[Alvarez,  2011],  and  as  such  is  an  important  area  of  study.  Three-dimensional,  Lagrangian  finite 
element  analysis  is  used  [Belytschko  et  al.,  2000].  One  long-term  objective  of  this  research  is  to 
have  the  capability  to  provide  uncertainty  quantification  and  robust  validation  for  an  intact  lower 
extremity  model  with  which  to  study  underbody  blast  effects.  For  this  reason  we  adopt  the  method 
of  Thacker  et  al.  [2007]  who  used  a  hierarchical  methodology  for  creating  a  numerical  model  of  the 
spine.  This  paper  documents  our  initial  efforts  in  building  the  computational  framework  required 
for  such  a  hierarchical  approach,  inspiring  both  smaller  and  larger  length  scale  studies  (i.e.  for  the 
tibia  and  full  leg,  respectively)  as  the  objects  of  future  work. 


2  Computational  Methodology 

Thacker  et  al.  [2007]  provide  a  useful  hierarchical  verification  and  validation  scheme  that  decom¬ 
poses  a  system  level  model  into  components,  subassemblies  and  assemblies.  Within  the  validation 
and  verification  scheme,  the  process  should  start  at  smaller  length  scales,  such  as  the  components, 
and  progress  to  the  last  stage  of  system-level  models.  This  approach  requires  simulations  and  ex¬ 
periments  to  be  conducted  at  each  level.  While  the  current  research  effort  adopts  the  approach 
proposed  by  Thacker  et  al.  [2007],  we  note  that  this  study  is  a  work- in-progress  and  is  still  in 
an  early  stage  of  development,  where  we  are  identifying  the  necessary  infrastructure  and  research 
teams  to  support  both  the  experimental  and  numerical  modeling  efforts  at  the  multiple  length  scales 
required.  Nevertheless,  we  begin  to  define  the  subdivisions  for  the  lower  extremities  for  this  ef¬ 
fort  and  to  also  help  pave  the  way  for  future  work.  The  lower  extremity  system  is  partitioned 
into  components  (ligaments,  muscles,  cortical  and  trabecular  bone,  skin),  subassemblies  (joints), 
assemblies  (motion  segments  consisting  of  muscles,  joints,  bones)  and  system  level  (whole  leg) 
models  as  depicted  in  Figure  la.  In  accordance  with  our  hierarchical  validation  and  verification 
strategy,  one  assembly  within  the  full  lower  extremity  model,  the  lower  leg  of  the  region  distal 
(below)  the  knee  is  considered  here. 

Twenty-eight  bones  are  explicitly  represented,  to  include  the  tibia,  fibula,  and  the  twenty-six 
bones  of  the  foot  [Hall,  2006].  The  tibia,  fibula,  talus,  and  calcaneus  are  constructed  as  layered 
thin-walled  structures  to  represent  cortical  and  trabecular  bone  types.  The  bones  of  the  mid-  and 
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(a) 


Velocity  (m/s) 


Figure  1:  (a)  A  three-dimensional  finite  element  volume  mesh  with  details  of  the  validation  hi¬ 
erarchy  for  the  lower  extremities  based  on  the  method  of  Thacker  et  al.  [2007].  The  system  is 
decomposed  into  (a)  components,  subassemblies,  assemblies  and  system-level  models,  (b)  Lower 
leg  model  geometry  in  assembly  model  of  leg  below  the  knee,  and  (c)  applied  velocity  boundary 
conditions  used  in  this  study. 


forefoot  are  not  segmented  into  different  bone  types  and  are  assigned  cortical  bone  properties 
throughout  and  have  the  ability  to  articulate.  Soft  tissue  and  a  separate  skin  layer  surround  the 
skeletal  structure.  The  soft  tissue,  hereby  referred  to  as  bulk  tissue  is  a  homogenization  of  mus¬ 
cles,  tendons,  connective  tissue,  blood  vessels,  etc.  using  an  approach  similar  to  that  proposed  by 
Cheung  et  al.  [2005].  The  bulk  tissue  also  fills  the  gaps  between  bones  so  that  joint  structures  such 
as  cartilage,  bursae,  synovial  fluid,  etc.  are  not  explicitly  represented.  Eight  ligaments,  however, 
are  explicitly  represented  in  the  model  and  are  detailed  in  Figure  lb.  The  proximal  (top)  end  of 
the  lower  leg  is  “potted”  in  a  ballast  mass  as  was  done  by  Bandak  et  al.  [2001].  The  bottom  of  the 
foot  contacts  a  layer  of  rubber  attached  to  a  steel  plate.  This  rubber-plate  arrangement  was  also 
used  by  Yoganandan  et  al.  [1995,  1996]  and  Bandak  et  al.  [2001]  as  a  means  of  applying  load  to 
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the  bottom  of  the  foot,  and  is  hereby  referred  to  as  the  impactor.  The  impactor  can  also  be  likened 
to  a  boot  sole  in  contact  with  the  floor  of  a  military  vehicle,  though  in  this  study  the  rubber  is  tied 
to  the  plate.  The  anatomical  geometry  for  the  model  was  obtained  from  Zygote  Media  Group,  Inc. 
[2012]  who  created  anatomical  geometry  from  CT  scans  of  human  specimens  who  anthropomet- 
ricly  fitting  within  the  50th  percentile  Zygote  Media  Group,  Inc.  [2012].  The  mesh  is  comprised 
of  1.53  million  10-noded  tetrahedral  elements  and  the  long  axis  of  the  tibia  is  aligned  with  the 
z-direction,  sometimes  referred  to  as  the  axial  direction.  Table  1  provides  details  on  the  material 
models  and  parameters  used. 

Both  cortical  and  trabecular  bone  failure  is  captured  within  the  simulations  using  simplistic 
failure  models  as  schematically  depicted  in  Figure  2.  In  tension,  the  material  behaves  as  a  linear 
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Figure  2:  Schematic  showing  the  elastic-fracture  constitutive  model  used  to  capture  fracture.  In 
tension,  first  there  is  an  linear  elastic  response  until  a  critical  stress,  o'..  is  reached  within  the 
element,  then  linear  stress  softening  occurs  until  the  element  critical  crack  opening  stretch,  £rc,  is 
reached.  In  compression,  element  fails  at  stress 


elastic  material  until  a  critical  stress,  o[:,  is  reached  within  an  element  at  strain  £e.  Then,  energy 
dissipation  associated  with  the  fracture  process  is  modeled  using  linear  stress  softening  until  the 
element  strain  reaches  a  critical  value,  £lc,  at  which  point  the  element’s  mass  is  also  removed  from 
the  system.  Note  that  the  critical  value  of  strain  includes  both  elastic  and  inelastic  contributions, 
i.e.,  £',  =  £e  +  £i,  where  £;  is  the  inelastic  dissipative  strain.  The  slope  of  the  stress  softening 
curve  should  be  scaled  so  that  the  energy  dissipated  by  the  degraded  element  equals  the  fracture 
energy  of  a  crack  passing  through  the  element,  thereby  rendering  the  solution  relatively  mesh  size 
independent  [Song  et  al.,  2008,  Donadon  et  al.,  2008,  Cervera  and  Chiumenti,  2006,  Iannucci  and 
Willows,  2006].  To  obtain  energy  equivalence  for  a  degraded  element  obeying  the  stress-strain  law 
shown  in  Figure  2,  the  following  relation  should  be  met: 

GfA  =  l-E£e£iVe  (1) 
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Table  1:  Compilation  of  various  constitutive  models  and  parameters  used  for  the  lower  leg  sim¬ 
ulation.  E  is  the  Young’s  Modulus,  v  is  Poisson’s  ratio,  K  is  the  bulk  modulus,  /i  is  the  shear 
modulus,  Kjc  is  the  mode-I  fracture  toughness,  o{.  and  ct?  are  the  critical  failure  stresses  in  tension 
and  compression,  respectively,  £,  is  the  inelastic  energy  dissipative  stretch,  p  is  the  density  and  m 
is  mass.  _ 


Anatomic  Component 

Material 

Model 

Material 

Properties 

References 

Cortical 

bone 

Rone 

Elastic  with 

maximum 

principal 

stress-based 

fracture  model 

E  =  19.0  GPa 

v  =  0.3 

Kjc  =  7  MPa  ^Jrn 
a[  =  132  MPa 
£,  =  0.022 
<JC£  =  -396  MPa 
p  =  1810  kg/m3 

Bandaketal.  [2001] 
Schuster  et  al.  [2000] 

Ural  et  al.  [2011] 

Currey  [2002] 

Trabecular 

bone 

Elastic  with 

maximum 

principal 

stress-based 

fracture  model 

E  =  300  MPa 

v  =  0.45 

Kic  =  0.4 

MPai/m 
&c=  1.5  MPa 
£,  =  0.230 
<7c£  =  -4.5  MPa 
p  =  600  kg/m3 

Schuster  et  al.  [2000] 

Cook  and  Zioupos  [2009] 

Ligaments 

Neo-Hookean 

K  =  250  MPa 
p  =  5.0  MPa 
p  =  1040  kg/m3 

Bandaketal.  [2001] 

Skin 

Neo-Hookean 

K  =  333  MPa 
p  =  6.7  MPa 
p  =  1040  kg/m3 

Bandaketal.  [2001] 

Bulk  Muscle 

Neo-Hookean 

K  =  333  MPa 
p  =  6.7  MPa 
p  =  1040  kg/m3 

Bandaketal.  [2001] 

Rubber  Layer 

Elastic 

E  =  1.0  GPa 

v  =  0.4 

p  =  1000  kg/m3 

Plate 

Elastic 

E  =  190.0  GPa 
v  =  0.28 
m  =  24.5  kg 

MatWeb  [2011] 

Bandaketal.  [2001] 

Knee  Ballast  Mass 

Elastic 

E  =  190.0  GPa 
v  =  0.28 
m  =  1 1.3  kg 

MatWeb  [2011] 

Bandaketal.  [2001] 

where  Gf  is  the  fracture  energy  of  the  material  per  unit  area,  A  is  the  fracture  area  of  a  crack 
passing  through  the  element,  E  is  Young’s  Modulus,  and  Ve  is  the  volume  of  the  element.  In  the 
current  study,  Gf  is  determined  using  the  elastic  fracture  mechanics  relation  Kjc  —  E'Gf,  where 
K/c  is  the  mode-I  fracture  toughness  of  the  material.  Therefore,  in  plane  stress  where  E'  =  E, 
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Gf  was  determined  to  be  2578  J/m2  and  533  J/m2  for  cortical  and  trabecular  bone,  respectively, 
using  the  values  listed  in  Table  1.  Both  A  and  Ve  in  Equation  1  are  mesh  parameters.  Assuming 
regular  tetrahedral  elements  in  which  all  edge  lengths,  a,  are  equal,  Ve  =  (vT/12 )a3.  Also  assuming 
a  flat  crack  face  passing  through  two  nodes  of  the  tetrahedron  and  bisecting  the  opposite  edge,  the 
fracture  surface  is  obtained  as  A  ~  (V2/4j<:r .  To  satisfy  Equation  1,  the  value  of  £,  is  adjusted  in  the 
elastic-fracture  material  model  for  both  cortical  and  trabecular  bone. 

In  the  current  mesh,  a  distribution  of  element  aspect  ratios  and  sizes  exist,  which  contributes 
to  a  range  of  element  volumes  Ve  (2.62  x  10-11  to  1.70  x  10-8m3  for  cortical  bone,  2.40  x  10~10 
to  9.36  x  10  8 in3  for  trabecular  bone)  and  corresponding  fracture  areas  A.  Given  this  variability, 
and  still  assuming  regular  tetrahedral  elements,  £,  has  a  range  of  0.193  to  0.022  for  the  smallest 
and  largest  elements  of  cortical  bone,  respectively.  For  trabecular  bone,  £,  has  a  range  of  1.683 
to  0.230  for  the  smallest  and  largest  elements,  respectively.  In  this  study  we  consider  how  these 
ranges  of  inelastic  strain  parameters  influence  the  predicted  values  of  stress  and  damage  in  the 
lower  extremities,  but  first  we  establish  a  benchmark  simulations  using  the  maximum  element 
length  scales.  The  inelastic  dissipative  strain  for  both  cortical  and  trabecular  bone  is  given  by: 

Cortical  bone: 


2 GfA  _  2(2578  J/m2) (9.72  x  10“6  m2) 
&cye  ~  (132  MPa) (1.70  x  l(T8m3) 


(2) 


Trabecular  bone: 

2(533  J/m2) (3.03  x  1(T5  m2)  _ 

F  —  -A _ A _ L  ~  n  230  131 

1  (1.5  MPa)  (9.36  xl0-8m3) 

Compressive  failure  in  both  cortical  and  trabecular  bone  is  captured  using  a  minimum  principal 
stress  failure  criterion  and  artificial  element  deletion:  Energy  is  not  dissipated  by  the  constitutive 
model  once  the  element  reaches  its  failure  criterion  and  its  mass  is  simply  removed  from  the  system 
upon  deletion.  As  shown  in  Figure  2,  no  stress  softening  curve  is  present  past  the  failure  criterion, 
<jcc,  because  the  elastic-fracture  material  model  used  in  this  study  currently  does  not  support  it. 
Adding  this  capability  will  be  the  focus  of  future  work.  occ  for  both  cortical  and  trabecular  bone  is 
determined  by  a‘:.  =  —  3  rr/. . 

A  variety  of  boundary  conditions  have  been  investigated.  Figures  lc  and  d  show  the  range  of 
velocity  profiles  applied  to  the  plate  in  this  study,  all  in  the  positive  z-dircction.  The  17.5  and  3.2 
m/s  velocity  profiles  (Figure  lc)  are  considered  baselines.  In  all  cases  the  plate  is  fixed  in  the  x 
and  y  directions.  A  gravity  setting  with  a  gravitational  constant  of  9.8  m/s2  has  been  applied  in  the 
negative  z-dircction.  The  ballast  mass  is  fixed  in  the  x  and  y  directions. 


3  Results  and  Discussion 

Figure  3  shows  the  evolution  of  minimum  principal  stress  predicted  by  the  computations  at  the  mid¬ 
tibia.  For  the  tibia  response,  both  tensile  and  compressive  stresses  were  examined  to  predict  failure; 
however,  the  tensile  component  is  small  in  magnitude  compared  to  the  compressive  component  and 
is  not  shown  in  the  figure.  The  compressive  stresses  reach  a  peak  value  of  approximately  -58  MPa 
for  the  17.5  m/s  case  (case  2)  at  approximately  3.2  ms;  for  the  3.2  m/s  case  (case  1),  peak  stress 
is  approximately  -38  MPa  at  3.8  ms.  The  percentage  of  bone  damage  is  also  shown  in  Figure  3. 
A  quick  estimate  of  damage  to  the  various  bones  can  be  obtained  from  the  percentage  of  mass 
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removed  from  the  system.  As  detailed  in  Section  2,  for  tensile  failure  the  mass  of  fully  degraded 
elements  is  zeroed,  while  for  compressive  failure  direct  element  deletion  is  used.  For  the  17.5  m/s 
case,  numerical  instabilities  occur  at  a  total  percentage  of  bone  damage  equal  to  11.5%;  for  the 
3.2  m/s  case,  numerical  instabilities  occur  at  a  total  percentage  of  bone  damage  equal  to  13.5%. 
This  total  bone  damage  will  be  decomposed  into  individual  components  and  discussed  in  more 
detail  shortly.  Interestingly,  for  case  2  the  foot  stays  in  direct  contact  with  the  rubber  surface  and 
experiences  significant  comminution  and  failure  -  the  decrease  in  stress  magnitude  in  the  mid  tibia, 
illustrated  in  Figure  3,  is  due  to  degradation  in  the  ability  to  carry  load.  However,  for  case  1,  while 
some  failure  is  predicted  much  of  the  stress  drop  in  the  tibia  is  due  to  unloading  of  the  lower  leg 
because  of  separation  of  the  foot-rubber  interface  (sometimes  referred  to  as  liftoff)  and  movement 
of  the  ballast  mass  at  the  top  of  the  lower  leg  assembly. 


Figure  3:  Minimum  principal  stresses  predicted  in  the  mid-tibia  for  the  3.2  m/s  and  17  m/s  loading 
cases.  Also  shown  is  the  mass  removed  from  simulations  as  a  result  of  failure. 


Figure  4a  and  b  show  the  bone  structures  at  the  end  of  the  simulations  for  the  two  load  cases 
shown  in  Figure  lc.  A  single  calcaneus  fracture  is  seen  for  case  1,  while  for  case  2,  massive 
comminution  of  the  calcaneus  is  observed.  Fracture  of  the  talus,  distal  tibia  and  comminuted 
fracture  of  the  calcaneus  is  also  predicted.  In  both  cases,  no  mid-shaft  tibia  fractures  are  predicted. 
Interestingly,  the  arrows  in  Figures  4a  and  b  shows  an  example  of  how  the  dynamic  fracture  process 
is  different  between  the  two  different  loading  rates.  In  case  1,  at  the  lower  strain  rate,  no  dynamic 
crack  branching  was  predicted,  however  in  case  2  crack  bifurcation  is  predicted.  This  type  of 
information  would  be  useful  to  the  armor  designer  or  medical  team  because  these  injuries  require 
different  treatment  approaches.  In  case  1  a  relatively  simple  fracture  is  predicted,  while  at  higher 
rates  bone  fragmentation  is  expected. 

Figure  5a  and  b  shows  the  global  minimum  stress  in  the  axial  direction  in  the  tibia,  talus  and 
calcaneus  as  a  function  of  time.  For  case  1  the  compressive  stress  in  the  tibia,  talus  and  calcaneus 
reach  approximately  -60  MPa  within  approximately  2.0  ms  of  impact.  At  approximately  2.0  ms 
the  rate  of  growth  of  the  compressive  stress  within  the  talus  increases  and  stresses  reach  a  max¬ 
imum  value  of  approximately  -200  MPa  at  3.9  ms  followed  by  unloading.  The  increase  in  talus 
compression  stress  occurs  just  as  elements  in  the  talus  begin  to  undergo  deletion  due  to  material 
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Case  1 

3.2  m/s  in  2.6  ms 


Case  2 

17.5  m/s  in  6.1  ms 


Figure  4:  Predicted  fracture  morphologies  for  boundary  condition  (a)  3.2  m/s  and  (b)  17.5  m/s. 


failure  (at  2.0  ms).  The  tibia  response  also  shows  a  sharp  increase  in  axial  stress  level  at  approx¬ 
imately  2.6  ms,  followed  by  a  gradual  stress  increase  to  a  maximum  of  -375  MPa  at  4.5  ms,  and 
subsequent  unloading.  The  calcaneus  response  shows  a  dramatic  increase  in  compressive  stress  at 
approximately  5.2  ms  when  the  axial  stress  reaches  the  compressive  failure  criterion  of  -396  MPa. 
Typically,  these  conditions  indicate  that  a  crushing  mode  has  been  reached  in  the  bone.  However, 
given  that  the  foot  is  entering  into  an  unloading  state  (as  indicated  by  the  displacement  of  the  bal¬ 
last,  the  dashed  line  in  Figure  5a)  numerical  instabilities  might  play  a  role  in  these  results,  which 
should  be  therefore  regarded  with  some  level  of  caution.  Despite  these  considerations,  for  both  the 
3.2  and  17.5  m/s  cases  damage  starts  first  in  the  talus.  Reflecting  waves  induce  localized  regions 
of  tensile  stress  where  cracks  can  initiate.  In  addition,  the  calcaneus  is  some  what  protected  by  its 
ability  to  rotate  counterclockwise  about  the  talus  in  the  sagittal  plane  while  the  tibia  is  able  to  with¬ 
stand  higher  levels  of  imposed  distal  displacements  as  the  strain  is  distributed  over  its  considerable 
length.  Conversely,  the  talus  is  tightly  constrained  between  the  rotating  calcaneus  and  the  tibia. 
For  the  17.5  m/s  case,  the  response  at  times  less  than  1.5  ms  is  qualitatively  similar  to  the  lower 
rate  case,  with  analogous  patterns  of  stress  growth  in  the  tibia,  talus  and  calcaneus.  Unfortunately, 
numerical  instabilities  are  observed  in  all  three  regions  as  the  transient  analysis. 

v  =  3.2  m/s  v  =  1 7.5  m/s 


2.5 


o> 

15  I 


N 

0.5 


Figure  5:  Maximum  compressive  stress  in  the  axial  direction  as  a  function  of  time  for  case  of  (a) 
3.2  and  (b)  17.5  m/s. 


Figure  6a  and  b  shows  the  principal  stresses  in  the  plantar  surface  of  the  calcaneus  and  inferior 
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talus  for  case  1  and  case  2,  respectively.  This  is  important  to  examine,  considering  the  large 
amount  of  failure  observed  in  these  regions.  One  interesting  aspect  of  the  response  is  the  presence 
of  tension-compression  components  resulting  from  the  combined  loading  in  these  regions.  In  the 
inferior  talus  the  rate  at  which  compressive  stresses  grow  is  faster  then  the  growth  rate  for  tensile 
stresses  up  to,  approximately  2  ms  into  the  transient  when  fracture  initiates  in  the  trabecular  region 
of  talus.  At  2.6  ms  the  acceleration  of  the  applied  velocity  boundary  conditions  reverses,  thus 
displacements  are  still  increasing  in  the  axial  direction,  but  the  rate  of  displacement  of  the  floor 
plate  decreases  with  time.  Case  2  shown  in  Figure  6b  shows  similarities  in  the  loading  of  the 
inferior  talus  up  to  the  time  of  initiation  of  fracture.  It  is  important  to  keep  in  mind  that  at  this 
time,  rate  dependent  material  properties  of  the  bone  deformation  and  fracture  response  are  not 
incorporated  into  the  model  so  any  stiffening  of  the  model  response  derives  from  inertial  effects 
alone.  By  approximately  2.4  ms,  fracture  morphology  shows  severe  comminution  (Figure  4).  An 
interesting  aspect  of  the  response  is  the  rate  at  which  fracture  occurs  at  a  given  location.  For 
example,  as  seen  in  Figure  6b,  at  approximately  3.6  ms  the  principal  stress  in  the  inferior  talus 
rapidly  grows  until  the  critical  stress,  Glc,  is  reached  and  the  failure  process  begins  to  dissipate 
energy.  The  curves  in  Figures  6a  and  b  represent  data  extracted  from  one  particular  region  of  the 
respective  model  components. 


v  =  3.2  m/s  v  =  17.5  m/s 


Figure  6:  Principal  stresses  in  the  plantar  surface  of  the  calcaneus  and  inferior  talus  for  both  (a) 
case  1  and  (b)  case  2. 


Figure  7  shows  the  damage  predicted  for  the  17.5  m/s  impact  case  at  subsequent  instants  of 
the  simulation  transient.  The  color  contours  indicate  if  the  damage  is  associated  with  tension 
or  compression,  and  whether  it  is  predicted  to  occur  in  the  trabecular  or  cortical  bone.  Damage 
initiates  in  the  talus  and  begins  at  approximately  1.45  ms.  At  2.50  ms,  a  mix  of  trabecular  and 
cortical  bone  failure  is  observed.  For  both  loading  conditions,  mid-tibia  failure  is  not  observed, 
suggesting  that  direct  axial  loading  for  the  lower  extremity  subjected  to  underbody  blast  mainly 
affects  the  ankle  complex.  However,  the  extrinsic  effects  such  as  lateral  seat  loading  to  the  mid 
lower  leg  giving  rise  to  off-axis  loading  may  significantly  alter  the  predicted  fracture  morphologies. 

Once  elements  are  fully  damaged,  they  are  removed  from  the  simulation  using  the  element 
deletion  method.  While  primitive  in  approach,  element  deletion  helps  to  keep  the  simulation  nu¬ 
merically  stable  and  also  provides  a  simple  and  quick  means  to  estimate  damage.  Figure  8  shows 
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v  =  17.5  m/s 


t  =  1.45  ms 


t  =2.50  ms 


t=3. 71  ms 


Gold  =  tensile  trabecular  failure 
Teal  =  compressive  trabecular  failure 
Red  =  tensile  cortical  failure 
Purple  =  compressive  cortical  failure 


Figure  7:  Damage  predicted  for  the  17.5  m/s  impact  case  over  various  times 


the  percentage  of  mass  removed  from  the  simulation  for  the  talus,  which  is  interesting  to  examine 
because  this  is  where  damage  is  first  observed.  This  is  not  the  most  accurate  means  to  evaluate 
damage,  because  this  measure  only  accounts  for  elements  fully  damaged.  In  fact,  many  elements 
using  the  elastic-fracture  model  are  in  the  “transition”  phase  of  failure,  dissipating  energy  associ¬ 
ated  with  the  unloading  portion  of  the  curve  shown  in  Figure  2.  For  the  3.2  and  17.5  m/s  cases, 
although  the  time  scales  of  damage  are  different,  similar  amounts  of  mass  are  removed  from  both 
tensile  and  compression  failure,  i.e.,  13  and  31.5  percent,  respectively.  However,  for  the  17.5  m/s 
case  cortical  failure  is  observed,  whereas  for  the  3.2  m/s  only  trabecular  failure  in  the  talus  was 
observed.  This  reveals  an  interesting  aspect  of  lower  extremity  failure  mechanisms  directly  rele¬ 
vant  to  the  investigation  of  the  injury  process  for  sub-critical  loading  that  cause  no  cortical  failure. 
The  cortical  bone  has  higher  strength  and  is  thought  to  directly  correlate  with  the  ability  of  the 
lower  extremities  to  carry  load.  The  lower  strength  of  the  trabecular  bone  makes  it  vulnerable  to 
dynamic,  non-equilibrium  stresses  that  can  activate  cracks.  This  is  an  important  deviation  from 
traditional  biomechanics  that  tend  to  study  the  bone  system  in  equilibrium  when  the  response  of 
the  cortical  shell  is  primarily  responsible  to  determine  the  load  bearing  capacity  of  the  structure. 


4  Effect  of  Critical  Crack  Stretch  Parameter 

Since  the  approach  used  in  this  study  is  one  that  is  mesh  dependent,  an  initial  parametric  study 
was  conducted  to  examine  the  effect  of  parameters  used  within  the  elastic  fracture  model.  Figure 
9  shows  the  maximum  compressive  principal  stress  predicted  by  the  model  plotted  in  terms  of  the 
tensile  critical  inelastic  stretch,  £,  for  the  17.5  m/s  loading  condition.  Also  shown  is  the  corre¬ 
sponding  total  mass  removed  from  the  system  as  a  result.  The  predicted  maximum  compressive 
principal  stress  is  diminished  for  low  values  of  critical  inelastic  stretch.  The  smaller  the  critical 
stretch,  the  less  energy  is  dissipated  during  the  failure  process,  leading  to  more  elements  becoming 
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v  =  3.2  m/s  V  =  17.5  m/s 


Figure  8:  Percentage  of  mass  removed  from  the  simulation  for  talus  for  the  (a)  3.2  m/s  and  (b)  17.5 
m/s  cases.  The  talus  is  interesting  to  examine  because  this  is  where  damage  is  first  observed. 


fully  damaged,  or  more  total  mass  removed.  A  significant  increase  in  maximum  compressive  prin¬ 
cipal  stress  occurs  when  the  critical  inelastic  stretch  is  increased  up  to  0.05.  Above  this  value  of 
stretch,  there  seems  to  be  a  less  dramatic  effect  on  the  principal  stress  and  a  lower  effect  on  damage 
reduction.  In  terms  of  fracture  surfaces,  a  lower  fracture  area  is  created  with  increasing  inelastic 
stretch.  In  the  future,  we  hope  to  modify  the  constitutive  law  to  take  into  account  the  distribution 
of  element  length  scales.  This  development  will  allow  the  model  to  account  for  a  distribution  of 
critical  inelastic  stretch  values,  thus  reducing  mesh  dependency  effects. 


5  Conclusions  and  Future  Work 

This  paper  has  described  a  research  effort  to  develop  a  hierarchical  modeling  approach  for  the 
lower  extremities  subjected  to  military  underbody  blasts  and  has  specifically  focused  on  devel¬ 
oping  a  finite  element  model  of  the  lower  extremities  undergoing  high  strain  rate  blast-induced 
deformation.  Computations  highlight  the  importance  of  understanding  the  initiation,  propagation 
and  coalescence  of  bone  fracture.  At  the  assembly  level,  a  lower  leg  consisting  of  many  different 
tissues  was  developed  and  evaluated.  Similar  to  past  experimental  research,  we  numerically  ob¬ 
serve  calcaneal,  talar,  midfoot  fractures,  as  well  as  ankle  rotation  before  any  tibia  or  fibula  fracture 
is  observed.  This  study  has  some  considerable  limitations  including  the  simplistic  approach  used 
to  study  the  fracture  process,  the  material  properties  that  neglected  rate  dependence  effects,  and  a 
number  of  anatomic  details  that  were  not  represented  but  may  be  important.  In  future  work,  com¬ 
ponents  within  the  assembly,  such  as  the  tibia  will  be  refined  by  using  sub-scale  models  of  trabec¬ 
ular  bone  at  the  micro  structural  level  and  will  capture  micro-fracture  mechanisms  at  the  trabecular 
level.  For  that  purpose,  we  will  make  use  of  the  Discontinuous  Galerkin  finite  element  method 
for  the  computational  model  [Radovitzky  et  al.,  2011].  Furthermore,  a  visco-elastic,  visco-plastic 
formulation  for  the  constitutive  response  of  cortical  and  trabecular  bone  [Johnson  et  al.,  2010]  will 
be  implemented,  to  allow  for  rate  effects  in  the  deformation  response  of  the  tissue.  To  better  un- 
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Figure  9:  Maximum  compressive  principal  stress  plotted  in  terms  of  varying  tensile  critical  in¬ 
elastic  stretch,  e,  for  the  17.5  m/s  loading  condition.  Also  shown  is  the  corresponding  total  mass 
removed  from  the  system  as  a  result. 


derstand  how  the  open  cellular  structure  of  trabecular  bone  affects  its  macroscopic  response  and 
failure  behavior,  a  meso-scale  model  will  be  developed.  This  hierarchical  approach  will  be  further 
extended  to  other  model  components  to  obtain  a  more  accurate  and  biofidelic  representation  of  the 
lower  leg  assembly. 
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